Prompt neutrino fluxes from atmospheric charm 
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We calculate the prompt neutrino flux from atmospheric charm production by cosmic rays, using 
the dipole picture in a perturbative QCD framework, which incorporates the parton saturation 
effects present at high energies. We compare our results with the next-to-leading order perturbative 
QCD result and find that saturation effects are large for neutrino energies above 10 6 GeV, leading 
to a substantial suppression of the prompt neutrino flux. We comment on the range of prompt 
neutrino fluxes due to theoretical uncertainties. 
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I. INTRODUCTION 

Atmospheric neutrinos are produced in interactions of 
cosmic rays with Earth's atmosphere. The observation 
of low energy [E v ~ GeV) atmospheric neutrinos, their 
flavor-dependent interactions, and their path length de- 
pendence [l], Q has confirmed the existence of neutrino 
flavor transformation, and therefore the most fundamen- 
tal property of the neutrinos: that they are not massless. 
These observations have provided a remarkable source of 
information about mass and mixing parameters of neu- 
trinos. 

Atmospheric neutrinos are also a background to other 
sources of neutrinos, such as cosmogenic neutrinos pro- 
duced in interactions of cosmic rays with the background 
radiation Q and directly from sources such as active 
galactic nuclei and gamma ray bursts [|| • Observation of 
neutrinos coming from these distant sources would pro- 
vide valuable information about the particle production 
mechanism in astrophysical sources. 

Neutrino interactions in the Earth and in the atmo- 
sphere could also serve as unique probes of physics be- 
yond the Standard Model [f| . It has been recently sug- 
gested that atmospheric neutrinos, in their interactions in 
the Earth, could produce supersymmetric particles if the 
neutrino energies are sufficiently high. In Ref. [6j], Ando 
et al. have suggested that the high energy atmospheric 
neutrino flux may be large enough to produce quasi- 
stable charged particles that are potentially detectable 
in the IceCube Q neutrino detector. As a background 
to high energy sources or as a flux to produce exotic par- 
ticles in the Earth, it is useful to re-evaluate the high 
energy component of the atmospheric neutrino flux. 

The atmospheric fluxes of neutrinos at low energies 
have been extensively studied [1, [t| Qil [H| ■ They arise 
mainly from the products of charged pion and kaon de- 
cays. As energies increase, the decay lengths of the 
mesons become longer than their path lengths in the 
atmosphere [l2j . suppressing the production of neutri- 
nos. Other, shorter lived hadrons are also produced 
at high energies. They too contribute to the neutrino 
flux, especially from the "prompt" decay of charmed 
mesons. The energy dependence of these prompt neu- 
trinos is less steep than the "conventional" neutrino flux 
from pion and kaon decays. The energy at which the 



prompt neutrinos become the dominant atmospheric neu- 
trino component depends on the details of the mech- 
anism for charm production in proton-air collisions at 
high energies. Charm contributions to the atmospheric 
lcpton fluxes have been evaluated analytically and semi- 
analytically QJ, Q M, EE E3 • 

There are renewed efforts 
include charm production with the dual parton model in 
Monte Carlo simulations of air showers as well [HI, E|| . 

For vertical neutrino fluxes, the cross-over between 
conventional and prompt dominated fluxes occurs in the 
energy ran ge o f 10 5 -10 6 GeV for the calculations of Refs. 
[li], [Til, 17 1, and the cross-over energy increases with 
zenith angle. For energies above ~ 10 6 GeV, the domi- 
nant contribution to charm production comes from glu- 
ons, where saturation effects [2(| due to dense, interact- 
ing gluons in the nucleus become important. We evaluate 
the prompt neutrino flux using perturbative QCD in the 
dipole framework, taking these effects into account. We 
study the theoretical uncertainties inherent in this ap- 
proach and compare with standard next-to-leading order 
perturbative QCD. The range of QCD-based predictions 
yields prompt neutrino fluxes that are unlikely to be large 
enough to produce a detectable number of exotic parti- 
cles of the type discussed in Ref. Q. 

We begin with a discussion of the cross section for 
charm production in the dipole picture in Section II. Us- 
ing the dipole picture results, we discuss the evaluation 
the prompt neutrino flux from charm decays in Section 
III. Our results and a comparison with the conventional 
fluxes are shown in Section IV. We also discuss uncer- 
tainties in the QCD approach, and compare our results 
with the uncertainty band of Ref. [6[ in Section IV. 



II. CROSS SECTION FOR CHARM 
PRODUCTION 



A. Charm production in perturbative QCD 

In the perturbative QCD approach, the dominant con- 
tribution to the charm cross section at high energies 
comes from the partonic subprocess gg — ► cc. The 
parton-level differential cross section for production of 
cc pairs in proton-proton collision, at the leading order 
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in the strong coupling constant, a s (/x 2 ), is given by 



d<JLO 
dxp 



dMl 



[x\ + x 2 )s 



(T gg ^ cS (s)G(xi,fi 2 )G(x2 1 ^ 2 ) (1) 



where x\ t 2 are the momentum fractions of the gluons, 
xp = x\ — X2 is the Feynman variable, G(x,[i 2 ) is the 
gluon distribution of the proton, and /i is the factoriza- 
tion scale. Given the charm-anticharm invariant mass 
M cg , the fractional momenta of the gluons, £1,2, can be 
expressed in terms of the the Feynman variable, Xf, 




± xf 



(2) 



Typically the factorization scale is taken to be of the 
order of 2m c . 

For the flux calculation we need the differential cross 
section as a function of incident proton energy (E p ) and 
final charm energy (E c ), convoluted with the incident 
cosmic ray proton flux. Clearly at high energies, given 
the relationship of Eq. ([2]), the charm cross section has 
dominant contribution when one gluon parton distribu- 
tion function (PDF) is at x\ ~ xp and the other gluon 
distribution is at 2; 2 <C 1. Since the gluon distribution 
cannot be measured directly, its value at very small x 
has large uncertainties, especially for the low factoriza- 
tion scale /1 ~ 2m c . The dipole picture gives a theo- 
retically motivated description of small x physics which 
can effectively take into account resummation of the 
large a s ln(l/a;) contributions [2l| to the evolution of the 
PDFs. Thus by using the dipole picture, we avoid the 
large uncertainty due to the unknown behavior of the 
gluon distribution at very small x. 

The dipole picture is most straightforwardly described 
in the DIS context, which we do next. We then elaborate 
how this is applied to hadron-hadron scattering. 



B. Dipole picture formalism in deep inelastic 
scattering 

In deep inelastic lepton-hadron scattering, the high Q 2 
virtual photon can penetrate the nucleon and probe the 
partonic degrees of freedom. This partonic interpreta- 
tion based on perturbative QCD is most relevant in the 
infinite momentum frame. The Q 2 -dependence of the 
nucleon structure function F 2 N (x 1 Q 2 ) is well accounted 
for by the DGLAP evolution equations [22j given some 
non-perturbative initial condition F 2 n (x,Qq). As noted 
above, at small x one needs to consider the resumma- 
tion of large logarithms In 1/x, which leads to the BFKL 
evolution equation [2ll ]. 

Another feature of the nucleon structure function 
in the DGLAP framework is the strong growth of the 
gluon density in the nucleon in the small x region. In 
the infinite momentum frame, because of the high nu- 
cleon and parton densities, quarks and gluons that be- 
long to different nucleons in the nucleus may recombine 



FIG. 1: The perturbative diagrams giving rise to scattering 
with a gluon of the 7* — ► qq fluctuation in deep inelastic 
scattering. 



and annihilate, leading to the recombination effect first 
proposed by Gribov, Levin and R ysk in (GLR) [2(| and 
later detailed by Mueller and Qiu [2a |. 

An alternative approach is to consider instead the in- 
teraction in the target rest frame (laboratory frame), 
where the virtual photon interacts with nucleons via its 
quark-antiquark pair (qq) color-singlet fluctuation [24[. 
If the coherence length of the virtual photon fluctuation 
is larger than the radius of the nucleus, l c > Ra, the qq 
configuration interacts coherently with all nucleons, with 
a cross section given by the color transparency mecha- 
nism for a pointlike color-singlet configuration [251 ] . That 
is, the cross section is proportional to the transverse sep- 
aration squared, r 2 , of the q and q. 

In the dipole picture, the cross section for the absorp- 
tion of a virtual photon in the small x region is dominated 
by the scattering of a gluon off the qq pair fluctuation of 
the virtual photon. The generic perturbative QCD di- 
agrams giving rise to the qq fluctuation are shown in 
Figure [I] The invariant mass of the incoming virtual 
photon-proton system at small x is related to the photon 
virtuality Q 2 by 



s = (q + p) 2 - 2p • q 



91 

x 



(3) 



where q and p are the four-momenta of the photon and 
the target nucleon, q 2 = —Q 2 and x = Q 2 /2p-q. Thus the 
region of small x corresponds to a high energy scattering 
process at fixed Q 2 . 

The imaginary part of the sum of the amplitudes in 
Figure [T] is related to the photoabsorption cross section, 
which has been calculated by Nikolaev and Zakharov [26| 
assuming that the size of the qq pair is frozen in the scat- 
tering process and that the one-gluon exchange process 
of Figure [1] dominates. The transverse cross section can 
be cast into an impact parameter representation 



a(j*N) 



J dz J d 2 r\t< T (z,r,Q 2 )\ 2 a q q N (x,r) , (4) 



where z is the Sudakov variable, defined to be the fraction 
of the qq pair momentum carried by the quark, and r is 
the variable conjugate to the transverse momentum of the 
quark, representing the transverse size of the pair. The 
function "$>t(z, r, Q 2 ) can be interpreted as the wave func- 
tion of the qq fluctuations of the virtual photon. Thus, 
\^t(z, r, Q 2 )\ 2 is the probability of finding a qq pair with 
a separation r and a fractional momentum z. It is given 
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FIG. 2: The perturbative diagrams giving rise to the scatter- 
ing of a gluon with the g — > qq pair fluctuation in hadronic 
collisions. 



for each quark flavor / with fractional charge e/ by [2f 



\H(z,r,Q 
2 a em N l 



2\|2 



(5) 



-r 



2tt 2 



[(z 2 + (l-z) 2 ) e 2 K 2 (er) +m 2 f K 2 (er)] 



-z)Q 2 



and Kn and Ki are modified 



where e 2 = z(l 
Bessel functions. 

The cross section for the high energy interaction of a 
small-size qq configuration with the nucleon, o"qqjv(r), can 
be calculated in leading-order perturbative QCD. In this 
approximation, one sets <j q qN(r) equal to [2 71 ] 



■pQCD 



— r 2 a s {p)xG(xi,^ 2 ). 



(6) 



This cross section is, as discussed above, proportional to 
the square of the size of the pointlike configuration as 
a consequence of color transparency in QCD. However, 
the singular behavior of the wave function and the strong 
scaling violation of the gluon distribution in the small-x 
region as r decreases can compensate the smallness of the 
cross section due to color transparency. 

Ultimately, gluon saturation effects need to be included 
for a more realistic cr q q^(r). One would then derive an 
approximate expression for the dipole cross section from 
theory, including saturation effects, and use experimen- 
tal data to determine incalculable parameters in this ex- 
pression. Before we turn to saturation and the types of 
functional forms used to fit the dipole cross section, in 
the next section we describe how heavy quark produc- 
tion in proton-proton scattering is treated in the dipole 
picture. 



C. Heavy quark production 

Heavy quark production in hadronic collisions can be 
obtained in the same formalism [2^, [2^, [3(3, [HJ . In this 
case, the dipole is produced from a gluon instead of a 
photon, so that the dipole can be in a color octet state. 
As shown in Figure^ there is now an additional diagram 
that contributes, in which the gluon interacts with the 
target before fluctuating to a dipole. 

The differential cross section for heavy quark produc- 
tion is [H 



dajpp -> QQX) 
dy 



Xl G{ Xl ^ 2 )a G ^ x {x 2 ^ 2 ,Q 2 



(7) 



where x\ and x 2 are the partonic momentum fractions, 
y = \ \rx[xxlx%) is the QQ pair rapidity and a G v~*QQ x is 
the partonic cross section calculated in the dipole model, 



(x,fi z ,Q') = J dzd z r\^(z,r)\'a dG (x,r) . 

(8) 

The probability of finding a QQ pair with a separation r 
and a fractional momentum z, is given by 



|*g(z,r,Q 2 = 0)| 2 = 
2ir 2 



(9) 



[{z 2 + (1 - zf) m 2 Q K 2 (m Q r) + m 2 Q K 2 {m Q r)] , 



where fi ~ 1/r is the factorization scale. For heavy quark 
production we have Q 2 — 0, so /i ~ uiq and e = tuq. 

The dipole cross section that describes the interaction 
of a heavy quark-antiquark pair from the fluctuation of 
a gluon with the target nucleon is given by [28[ 



T GQQ 
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{X2,r) = - [a d {x 2 , zr) + a d (x 2 , (1 - z)r) 



(10) 



where a d is the color singlet dipole cross section of Eq. 
(4). The first term corresponds to the quark-gluon 
(G—Q) separation zr, the antiquark-gluon (G—Q) sep- 
aration (1 — z)r and the quark-antiquark (Q — Q) sepa- 
ration r. This expression includes contributions from the 
three different color and spin states in which QQ can be 
produced [30| . 

Finally, to take threshold corrections for charm pro- 
duction at large x into account, the dipole cross section 
is multiplied with a factor (1 — X2) 7 [32|. We find this 
correction to be negligible for energies above 10 3 GeV. 



D. The dipole— proton cross section and saturation 

The dynamics of the scattering process at small x is, in 
principle, included in the dipole cross section. Thus, to 
compute the differential cross section da/dxF we must 
find the cross section for a cc dipole to scatter on the 
proton, including the effects of saturation. 

A simple model for saturation was proposed by Golec- 
Biernat and Wiisthoff [Hj]. In their model, the dipole 
cross section is parametrized as 



gbw „ 1 
a d = (To 1 - e 



*0)/4 



where Q s is the saturation scale, 

Q s = Q s (x) = Q (x /x) x / 2 



(11) 



(12) 



with Qq = 1 GeV. The parameters A and xo in the 
above expressions were fitted to HERA data on the struc- 
ture function F2 and the diffractive structure function 
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This is a phenomenological model, constructed to give 
the right behavior of the dipole cross section in the two 
limits r — > and r — > oo. Eq. (jlip has a oc r 2 for 
small r, as implied by perturbative QCD, and a — > const 
for large r (this is the saturation property of the cross 
section), thus providing some insight into the physics of 
saturation. The simple parametrization also gave a good 
fit to the data, although it does not reproduce newer data 
as well [HI as it does the older data. 

One would like to calculate the dipole cross section rig- 
orously in perturbative QCD; however, it is not known 
how to fully include the effects of saturation. It is con- 
venient to study QCD evolution in Mueller's dipole for- 
mulation [35j ] . where the projectile contains a collection 
of color dipoles. It has been shown [36| that in the 
high-energy limit, the scattering process is equivalent to 
a stochastic reaction-diffusion process where there are 
fluctuations in the number of dipoles. These fluctuations 
may potentially have a large effect on the energy depen- 
dence of the amplitude and saturation scale. A full cal- 
culation should include these effects, but they were found 
to be small in the region of very small x [371 ] . In princi- 
ple one should also take into account the complicated dy- 
namics of the color glass condensate [H, HH, Ho, EH . This 
is described by the functional intcgro-diffcrcntial Jalilian- 
Marian, Iancu, McLerran, Weigert, Leonidov and Kovner 
(JIMWLK) equations [39|], or equivalently by Balitsky's 
infinite hierarchy of coupled differential equations for the 
expectation values of Wilson lines [I(| • 

A much simpler equation which includes saturation 
was obtained by Balitsky [Z(| and Kovchegov [13] in the 
particular case where the target is a large nucleus. This 
equation is known as the Balitsky-Kovchegov (BK) equa- 
tion and although it can be derived within the dipole 
framework, it turns out to represent a specific mean-field 
approximation to the Balitsky- JIMWLK equations. The 
BK equation is, like the BFKL equation, a leading loga- 
rithmic evolution equation in ln(l/a;). The BFKL equa- 
tion, however, is a linear equation, while the BK equa- 
tion is similar in structure to the GLR [2(| and Muller- 
Qiu [HI equations, and can be written as the BFKL 
equation modified by a non-linear term. This reduces 
the power growth of the gluon distribution as x — > 0, 
which has been established by both numerical and ap- 
proximate analytical studies, see e.g. Ref. [37j and refer- 
ences therein. 

The dipole cross section is obtained from the solution 
of the BK equation, which can be solved numerically. We 
will instead use an approximate result [43| , which consists 
of a matching of approximate analytic solutions of the BK 
equation in the two regions of dipole size r 1/Q S and 
r <C 1/Qs) where the equation simplifies. 

In the large r region the solution approaches a fixed 
saturated value as r — ► 00 j44|. For r <C l/Q s the effects 
of the non-linearity in the BK equation are small, and the 
equation reduces to the BFKL equation; the solution is a 
saddle point solution of the BFKL equation, subject to a 
saturation condition 431 . The two solutions are matched 



Ref. Mo 7 a 



A 



•f.'o 



(Jo (mb) 



[43] 0.7 0.627 0.253 0.267 x 10~ 4 25.8 
[45] 0.7 0.627 0.175 0.19 x 10" 6 37.3 
[46] 0.7 0.738 0.220 0.163 x 10~ 4 27.3 

TABLE I: Parameter values in the dipole cross section for- 
mulas in Eqs. |l3j) and 111). In Refs. and [11, 7 3 is 
calculated, while in Ref. [461 ] . it is a fit parameter. 



at an intermediate scale rQ s = 2. The resulting model is 
what we will refer to as the dipole model (DM). 
The dipole cross section is given by 



a d (x,r) = a J\f(rQ s ,Y), 



(13) 



where 00 is a constant, and Af is the forward dipole scat- 
tering amplitude obtained from the BFKL or BK equa- 
tion I43I1. 



Af(rQ s ,Y) 



27eff(x,r) 



for r < 2 



Ml (I) ■ "" 1 ■ - 1 1.4) 

1 - cxp [-a\n 2 (br)] , for r > 2 . 



Here r = rQ s , Y = ln(l/x) is the rapidity, and again the 
saturation scale is defined in Eq. (fT2")) with Qq = 1 GeV. 
Furthermore, 



7eff(a;, r) = 7 S + 



ln(2/r) 
k\Y 



(15) 



is the "effective anomalous dimension," and 7 S and n are 
theoretical parameters calculated from the BFKL equa- 
tion, with numerical values 7 S = 0.627 and k = 9.94. 
Note that this is a perturbative QCD result and not an 
ad hoc model, although it is obtained by an approximate 
solution of the BK equation. 

The free parameters in the model are Afo, <Jq, A and 
xq . In Ref. [43j , the first of these was chosen to take the 
value 0.7. The exponent A specifies the power behavior 
of the saturation scale with x, and xq is the value of x 
where the saturation scale is 1 GeV. Furthermore, a and 
b are matching coefficients to be chosen such that the 
dipole amplitude and its derivative with respect to r are 
continuous at r = 2. Wc find 



ln(l-AAo) 



ln 2 (l-A/" )^ 



b=-(l-JV )- 



(16) 
(17) 



Note that the amplitude is a function of r and x only in 
the combination indicated, r = rQ s {x), except for the ge- 
ometric scaling breaking term in the effective anomalous 
dimension which contains the rapidity. 

The fitted parameter values from three different fits 
to HERA data on the deep inelastic structure function 
F 2 at small x [H, [H, Si] are shown in Table Q] Note 
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that in all cases Ao was fixed at No = 0.7. The first row 
shows the original parameter values obtained in Ref. [43[ . 
This was a three-flavor fit and is therefore not suitable 
for our calculation, but it has been extended to include 
charm [45j . giving the values in the second row. Finally, 
the third row shows a more recent fit by Soyez [46[ , which 
also includes charm. In this fit the parameter 7 S was 
taken as a free parameter, which gave a better fit to the 
data, with a larger value of "f s and a smaller value of 
A. This is quite interesting since a reduction of these 
parameters is exactly what is expected when including 
higher order logarithmic corrections to the BFKL kernel 
in the BK equation [47[ . 

These models take into account only the leading expo- 
nential ^-dependence of the saturation scale, and there 
are lar ge s ub-asymptotic corrections to the energy depen- 
dence 1481. 



]nQ 2jy ) = ^xM Y _J_ ]llY 



3 / 2vr 1 

tFV ax" (is) VY 



+ OQ./Y), (18) 



where x(t) is th° BFKL characteristic function and 7 S = 
0.627. The models discussed in this section keep only the 
leading term in this expression. Using the full expression 
could potentially change the energy dependence of the 
cross section substantially, but to incorporate this in this 
dipole model would require introducing more parameters 
and performing a new fit to all the data. 

In the dipole model results that follow, the DM results 
shown use the parameters of Soyez [46J shown in Table U 
and the paramctrization of equations (|13H14p . 



E. Nuclear effects 

In a dipole framework there are two possible ways to in- 
clude nuclear effects suggested in the literature: modifica- 
tion of the saturation scale, e.g. as pro posed by Armesto, 
Salgado and Wiedemann (ASW) [4!| (see also I50jl for 
another approach), and the Glauber-Gribov [5l|,[52] for- 
malism. In the former case, the nuclear effects are ac- 
counted for by geometric scaling, simply scaling the sat- 
uration scale for a nucleus A according to 



Ql,A 



Q s.j. 



AirRl 



\ nR 



(19) 



1.12AV3 



where R p is the proton radius, R A 
0.86^4 -1 / 3 fm is the nuclear radius, and 5 is a free param- 
eter to be fitted to data. ASW find S = 0.79 by fitting 
to 7*^4 data at small x. The proton radius is related to 
ao in the dipole cross section through cr — 2nRp. 

In the Glauber-Gribov formalism, nuclear rescatter- 
ing is taken into account by integrating the dipole cross 



section for dipole-nucleus collisions over the impact pa- 
rameter, 



(20) 



<r$(x,r) = / d z ba2(x,r,b), 



where b is the impact parameter between the center of 
the dipole and the center of the nucleus. The expression 
for the 6-dependent cross section is given by 



aj(x,r, b) 



1 — exp 



1 



AT A (b)a p d ( Xl r) 



(21) 



where a p d (x,r) is the dipole-proton cross section given 
in Eqs. (fT3")) and (|14p and T A (b) is the nuclear profile 
function, 



T A {b) = / dzp A {z,b), 



(22) 



where p A is the nuclear density, and T A is normalized so 
that 



J d 2 bT A (b) = 1. 



(23) 



This model has e.g. been used in Ref. [53[ with a Fermi 
distribution for p to compute nuclear structure functions 
with good results. 

We compared the Glauber-Gribov model with a Gaus- 
sian distribution for p to the ASW method and found 
that for the relatively light air nuclei, these two methods 
give very similar results (within 10%). We will in the 
following use the simpler ASW method. 

Predictions from the framework described above have 
been tested against data. For deep inelastic structure 
functions this was done in Refs. [43|, |46| . The ratio of 
DIS on nuclei to DIS on deuterons was calculated and 
compared to E665 data in Ref. [54J, and total cross sec- 
tions for TP, 7^4, pp, and pA were calculated in [3l[ using 
the parameters from [45| (second row in Table and 
the jp and pp results were compared to data with good 
agreement. There have been no tests in the energy range 
probed by cosmic rays. However, the LHC will begin to 
access these energy scales shortly. 



F. Fragmentation of charm quarks 

Our earlier analytical calculation [l4[ did not take frag- 
mentation of the charm quarks into charmed hadrons into 
account, but simply took the hadron to have the same en- 
ergy as the charm quark. In Ref. [l6| , fragmentation was 
taken into account by decreasing the momentum fraction 
of the hadron to an average lower value. In this paper we 
take fragmentation into account by using fragmentation 
functions. 

For our comparison without fragmentation, we use up- 
dated hadron fractions [55| 

f D0 = 0.565, f D+ = 0.246, f D + = 0.080, / Ac - 0.094 



(24) 
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Hadron h 


N h 


£/, 


D° 


0.694 


0.101 


D+ 


0.282 


0.104 


Dt 


0.050 


0.032 


A+ 


0.00677 


0.00418 



TABLE II: Parameters in the LO Kniehl and Kramer frag- 
mentation model 15711. 



these functions give the fragmentation fractions in the 
KK model, and these are quite different from the val- 
ues cited above: for the LO fit we obtain f D a = 0.745, 
f D + = 0.296, f D + = 0.125, and / Ac = 0.063. KK also 
perform a NLO fit. The NLO values decreases the nor- 
malization of the calculated flux by about 10%, but as 
our calculation is a LO calculation it is more consistent 
to use the LO fit. These fragmentation fractions do not 
add to one, an indication of one of the theoretical uncer- 
tainties. 



where fh is the fraction of fragmentation of c — > h. These 
newer values are somewhat different from the values used 
in [3, EH, EH] ; this increases the computed flux by about 
20%. 

In general the cross section for hadron production in- 
cluding fragmentation is obtained from the cross section 
for charm production as 



da(pp —* hX) 
dEh 



dE c da(pp -> cX) , 



dE c , 



D n c (Eh/E c ), 



(25) 

where D^{z) is the fragmentation function for c — > h. 
This can be written in terms of momentum fractions as 



da{pp —* hX) f 1 dz da(pp — > cX) 



dxi 



dx r 



D h c {z), (26) 



where z — Eh/E c , x c — E c /E p , and xe = Eh/E p . At 
high energy the momentum fraction xe — %F, or, for the 
charm cross section, x c ~ Xp- 

We use both the older Peterson fragmentation func- 
tion [56j and the recent parametrization by Kniehl and 
Kramer (KK) in Ref. [57j ■ The Peterson function is given 
by 



(27) 



where e = 0.05 is a fitted parameter [5 51 ] - common for 
all mesons, and Nh = fhN is a normalization constant 
where N is given by the condition 



E 

h 



dzD^(z) = 1, 



(28) 



assuming that the shape of is independent of the 
hadron h, and the fragmentation fractions are given 
in Eq. (|24|) . The calculation without fragmentation 
amounts to taking fragmentation functions D^(z) = 
hS(l-z). 

The KK fragmentation function has the form 



h x(l-x) 2 
D h c = N h - 



[(l-x) 2 + e h xY 



(29) 



with the parameters given in Ref. [57| , which we show in 
Table HH 

The Kniehl-Kramer fragmentation functions have nor- 
malization factors fitted to the data. The integrals of 



G. Theoretical uncertainties in the charm pair 
cross section 



Because the charm quark mass is of order 1 GeV, there 
are in principle large uncertainties in the charm pair pro- 
duction cross section |58j . In perturbation theory us- 
ing parton distribution functions, the charm cross section 
predictions can vary by more than an order of magnitude 
depending on the charm quark mass, number of flavors 
and choice of scales and PDFs. The dipole approach, 
with the fit to DIS data then translated to hadron scat- 
tering, mitigates the uncertainty. Beyond the total cross 
section, one is interested in the energy distribution of the 
charmed quark. 

To investigate the sensitivity of the charm differential 
cross section to the choice of parameters, we vary them 
as follows: We use the parameters of Ref. [46| for the 
dipole cross section (the fit of Ref. (4f| gives very similar 
results). We vary the PDF by taking the MRST 2001 
LO [I| or the CTEQ 6L gluon distributions and 
we vary the factorization scale between fip — 2m c or 
fip = m c , where the charm quark mass is varied between 
m c = 1.3 GeV and m c — 1.5 GeV. In each of the listed 
cases the former choice is what we use as our "standard" 
curves below. In Figure [3] we show a representative set 
of predictions for the differential cross section da(pA — > 
cc)/ dxp for A — 14.5, the average nucleon number of air, 
and an incident proton energy of 10 9 GeV. The parameter 
combinations that are not shown in the plot give results 
that fall between the upper and lower lines. 

We are also interested in the difference between the 
predictions of NLO QCD and the saturation prediction 
of the DM model. This is illustrated in Figured! where 
we show da(pA — > cc)/dxF at three energies using these 
two calculations. The NLO QCD cross section come from 
Ref. [3 (PRS). Note that the NLO QCD cross section 
increases with energy much faster than the DM cross 
section. For the lower energy E = 10 3 GeV the cross 
sections are comparable, but we shall see that because of 
the different energy dependence, the neutrino flux calcu- 
lated using NLO QCD is larger than the one calculated 
from the DM model. 
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The cascade consists of production and attenuation 
through interactions and decay of the particles. We fol- 
low the analytic approximation method used in Refs. 
[3 EE IH Isl for calculating the flux. In Ref. [jj] it 
was shown that this approximate solution agrees with a 
numerical Monte Carlo solution of the same equations. 

The flux is calculated as a function of the slant depth 
X, which is a measure of the amount of atmosphere tra- 
versed by the particle. It is defined as the integral of 
the atmospheric density along its path through the at- 
mosphere: 

/•OO 

X{t,ff) = / d£'p(h(£',6)), (31) 
Jt 



FIG. 3: Charm quark xf distribution in proton-air collisions 
at E p — 10 9 GeV, calculated in the dipole model described 
in the text, with the "standard" choices of m c — 1.3 GeV, 
factorization scale fi F = 2m c and the MRST2001 PDFs [g(}. 



where h(£,8) is the height at distance from the ground I 
and zenith angle 9. A reasonable model for our purposes 
is an exponential atmosphere with [IE] 

p(h) = p exp(-h/h ), (32) 




with ho = 6.4 km and p = 2.03 x 1(T 3 g/cm 3 . The 
vertical depth of the atmosphere is X ~ 1300 g/cm 2 , 
while the horizontal depth is X ~ 36,000 g/cm 2 . We 
shall mostly be concerned with the vertical flux, 6 = 0, 
as the conventional flux is the smallest in the vertical 
direction. We will, however, show predictions for the 
flux in the horizontal direction as well. 

The general form of the cascade equation for the flux 
<t>j = <j>j (X, E) of particle species j at energy E and slant 
depth X is 



dX ~ \- Xdcc+Z^ fc ^b 



(33) 



FIG. 4: The NLO QCD pA -> ccX differential cross sec- 
tion as a function of Feynman xf for PRS [l4| compared to 
the dipole model (DM) result for incident proton energies of 
10 3 , 10 6 , 10 9 GeV. The thicker lines are PRS and the thinner 
lines with the same color are DM at the same energy. 



where Xj is the interaction length, A^ ec is the 



length, and S(k 

by 



S(k -> j) 



j) is the regeneration function. 



dE 



,(j> k {E') dn(k^j;E',E) 



X k (E') 



dE 



decay 
given 



(34) 



III. CALCULATION OF NEUTRINO FLUXES 

The lepton flux at sea level is calculated by solving 
the coupled set of differential equations that describes 
the cascade in the atmosphere initiated by the incident 
cosmic ray nucleons. We use the primary nucleon flux 
parametrization with a knee from Ref. [15| | : 



4>n(e) 



1.7 E 

174 E 



-2.7 



-3 



for E < 5 ■ 10 6 GeV 
for E > 5 • 10 6 GeV, 



where the cosmic ray energy E is given in GeV and the 
flux (j) N {E) in cm- 2 s~ 1 sr- 1 (GeV /A)' 1 . 



For the case of production, 

dn(k -> j; E k , E : ) _ 1 da(kA -> jY, E kl E 3 ) 



dE, 



O-kA(Ek) 



dEi 



is the distribution of secondary hadrons and a k A is the 
total inelastic cross section for kA collisions. For the case 
of decays, 



dn(k -> j; E k ,Ej) _ 1 dT(k -> jY, Ej) 
dEj ~ fk dEj ' 



(36) 



The nucleon, meson, and lepton fluxes are described by 
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the equations 

d<t>N 
dX 
d<t>M 
dX 



X 



N 



+ S(NA -> NY) 



S(NA -> MY") - 



pdM(E) 



^ + 5(MA -> My) 

Am 



^ = 5>(m^) 



(37) 

(38) 
(39) 



M 



where I = fj,,f^,,y e and the mesons include unstable 
baryons: for prompt fluxes from charm M = D ± , D°, 
D°, Df,Af. In Eq. (ggj) d M = c/3jt is the decay length. 

The analytic solution relies on the approximate factor- 
ization of the fluxes into energy- and X-dependent parts. 
For the meson flux: 



HI 



dX 



pd 



M 



— — + Zmm-t— + Znm-t— (40) 
Am Am Ajv 



with 



->kh 



dE 



,(f> k (E',X,d) X k (E) dn(kA — > hY; E' , E) 
<t> k {E,X,9) X k {E') dE 



( 41 ) 

We now make the standard assumption that 
<f> k (E,X,9) = E~^<j) k {X,9), so that if the energy 
spectrum falls as _E~ 7_1 , we have 



-7-1 



X k (E) dnjkA -> hY;E',E) 

XkW) dE ■ 

(42) 

Eq. (|3T|) for the nucleon flux then has the solution 

<t> N (X,E) = 4,(E)e- x / A »( E \ (43) 

where <j){E) = 0(0, E) is the primary flux of nucleons on 
the atmosphere and Apf(E) is the nucleon attenuation 
length, defined as 



A N (E) 



1-Z NN (E)' 



(44) 



where Xn(E) is the interaction length of nucleons in the 
atmosphere. It is given by 



X N (E) = 



A 



N a pA (E)' 



(45) 



where A = 14.5 is the average atomic number of air, A^o 
is Avogadro's number, and u v a is the total nucleon-air 
cross section. We take the parametrization from [64( for 
this cross section, and the Monte Carlo result from (l5l | 
for Z NN (E). 

The meson fluxes are expressed in terms of the nu- 
cleon flux by solving the cascade equations separately at 
low and high energies, where the interaction and regen- 
eration terms and the decay terms, respectively, can be 



neglected. For the high energy flux we need the atten- 
uation lengths of charmed hadrons in the atmosphere, 
which we replace by the corresponding quantities for K- 
mesons. These are approximated by 



A M {E) 



A 



(Tpp(E) 



1 



N a pA (E) <JKp(E) 1 - Z KK {E) 



(46) 



As for nucleons, we take Zkk from J15J and u p a from 
[64j . The cross sections a pp and ok p are taken from [55[ . 

The final step is to obtain the lepton fluxes at high 
and low energies from Eq. (j39|) and the obtained meson 
fluxes, and interpolating between them for intermediate 
energy. This calculation is done in the limit X — > oo. The 
Z-moments for the three-body decay modes M — * IY 
are calculated using expressions in Refs. [H, [64|, and 
the lepton flux at intermediate energies is obtained by 
interpolating between the high- and low-energy solutions. 
In each of these regimes the meson fluxes are described 
by power laws 4>m{E) oc E~@ where (3 — 7 in the low- 
energy regime and /3 = 7 + 1 in the high energy regime, 
and 7 is the index of the primary nucleon flux. The 
higher power of energy in the high energy flux is due to 
the appearance of the gamma factor in the decay length 
in the denominator of the meson flux. 

The equations for the lepton fluxes then give 



- low _ 7 Z NM , (Tr s 
Pi - Z M l,~i+\Z 5 9n{E) 



(47) 



,hi g h ry z nm ln(A M /Ajv) £m , rj? s ,. Q % 
4>i =Z M l 1+ 2- ^ T—r-r ^(Pn{E), (48) 

where 6m, the critical energy for meson M, separates the 
low- and high-energy regions, where attenuation is dom- 
inated by decay and interaction. It depends on zenith 
angle, and is for the specific model of the atmosphere we 
use given by 



e M (9) 



mMC 2 h 

CT M 



(49) 



where h = 6.4 km is a scale parameter for the isother- 
mal height dependence of the atmospheric density [l5[. 
For relatively small angles, f(0) = 1/ cosd, but for an- 
gles near horizontal, the angular dependence is more 
complicated. To compute the horizontal flux, we fol- 
low the approach of Ref. [iH , leading to the replacement 
9 = 90° -> 9* = 84.45°. 

Further details of this procedure to solve the cas- 
cade equations semi-analytically are given, e.g., in Refs. 
[13 , [HI • Our treatment here adds the fragmentation of 
the charm quarks into charmed hadrons, meaning that 
we must compute separately the moments Z p h for each 
hadron M, including fragmentation functions in the cal- 
culation of the cross section. When fragmentation is ne- 
glected, we have the simple relation Z p h = fhZ pc , where 
fh is the fragmentation fraction. 
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log 10 E [GeV] 



FIG. 5: Prompt and conventional v^ + v^ fluxes in the vertical 
direction. The shaded band is the theoretical uncertainty 
band for the prompt flux calculated in this paper with the 
dipole model. The dashed line shows the conventional flux 
from Gaisser and Honda (GH) and the dotted line is the 
conventional flux calculated in Ref. [l5| ] (TIG). 



IV. RESULTS AND DISCUSSION 

Our result for the vertical muon neutrino plus antineu- 
trino flux from atmospheric charm is shown in Figure 
which shows the theoretical uncertainty band for the DM 
calculation, estimated as described above. For compar- 
ison the conventional neutrino fluxes [Til . [l5| from n- 
and -ftf-decays are also shown. We find that the verti- 
cal prompt muon neutrino flux becomes dominant over 
the conventional neutrino flux at energies between 10 5 
GeV and 10 5 5 GeV. 

The theoretical uncertainty due to choices of gluon 
distribution, charm quark mass, factorization scale, and 
other parameters in the dipole model result in the range 
of fluxes represented by the shaded area in Figure [5] The 
shape of the prompt neutrinos is only weakly dependent 
on the choice of parameters, but the overall normaliza- 
tion could vary by up to a factor of two in this model for 
charm production. 

We compare our result to three earlier calculations of 
the prompt neutrino flux: 

1. Thunman, Ingelman and Gondolo (TIG) [Hj]. This 
was the first perturbative QCD calculation and was 
done at the leading order (LO) in a s . It takes 
the fragmentation of charm quarks into account 
through Monte Carlo simulation using the Lund 
string mod el 1 651] implemented in the event genera- 
tor Pythia [66| . The small- a; PDFs are extrapolated 
with e.g., xG(x,fj, 2 ) ~ a;" 08 . 

2. Pasquali, Reno and Sarcevic (PRS) This re- 
sult uses the next-to- leading order (NLO) QCD re- 
sult of [67j with power law extrapolations of the 
small- a: PDFs. The PRS evaluation does not take 
fragmentation into account. We have therefore car- 




log 10 E [GeV] 



FIG. 6: Prompt muon neutrino fluxes obtained in perturba- 
tive QCD. The shaded area represents the theoretical uncer- 
tainty in the prompt neutrino flux evaluated in this paper, 
and the solid line in the band is our standard result. The 
dashed curve is the NLO perturbative QCD calculation of 
Ref. [3] (PRS), modified here to include fragmentation, the 
dotted curve is the saturation model result of Ref. il6j (MRS), 
and the dash-dotted curve is the LO perturbative QCD cal- 
culation of Ref. [II (TIG). 



ried out a simplified version of this calculation, tak- 
ing fragmentation into account in the same way 
as we did for the DM calculation: we compute 
the charmed hadron cross section in leading order 
QCD using KK fragmentation functions [57], and 
multiply with a K-factor K = a(NLO)/a(LO) w 
2. This reproduces the full NLO calculation of 
Ref. [IH at the parton level to an adequate accu- 
racy. 

3. Martin, Ryskin and Stasto (MRS) [ijj]. This calcu- 
lation takes fragmentation into account by assign- 
ing the neutrino a fixed fraction of the momentum 
of the mother meson, and is done using the sat- 
uration model of Golec-Biernat and Wiistoff [33| 
described above. 

We show the results from these other evaluations of the 
vertical muon neutrino plus antineutrino flux together 
with our uncertainty band in Figure [5] The theoretical 
uncertainties in the standard NLO QCD calculation of 
the charm cross section are the choice of the renormal- 
ization and factorization scales, the charm mass, and the 
small x behavior of the gluon distribution [58| . The im- 
pact of some of these uncertainties on the neutrino flux 
has been studied in Ref. [17| . 

The MRS curve in Fig. [6] is at the lower border of 
our DM uncertainty band. There is approximately a fac- 
tor of two between the MRS and the central DM results, 
coming from the different parameterizations of ad- The 
enhancement is also seen in calculations of photoproduc- 
tion of heavy quarks [5!| comparing the GBW model and 
the improved DM model of Eq. [TjJ] The DM cross sec- 
tion for charm pair production in pp collisions lies within 
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log 10 E [GeV] 

FIG. 7: The effect of fragmentation on predicted + 
fluxes. The solid lines are DM and NLO QCD (the latter 
multiplied by two to separate the lines) with Kniehl-Kramer 
fragmentation. The dashed lines are without fragmentation. 




Iog 10 £ [GeV] 

FIG. 8: Theoretical uncertainty band (blue, dark band), com- 
pared to the uncertainty band from Ref. [6j| (magenta, light 
band) with the overlapping region shown as the middle (light 
blue) band. The flux is scaled by E, in units of l/km 2 yrsr. 



the uncertainty band of Ref. [58| . 

The effect of quark fragmentation on the neutrino 
fluxes is rather large because fragmentation reduces the 
energy of the charmed hadrons. For a given hadron en- 
ergy, fragmentation effects require higher energy cosmic 
rays in the steeply falling cosmic ray flux. In Figure [7| we 
show the effect of including the KK fragmentation func- 
tions on both the NLO QCD and DM results. The NLO 
results are multiplied by a factor of two so that they can 
be distinguished easily from the DM results. The frag- 
mentation reduces the flux by between 60% and 70%, 
and thus it is an important effect to take into account. 
The Peterson fragmentation function results differ by ap- 
proximately 10% from the results shown in Figure [3 We 
include the uncertainty in fragmentation model in our 
estimate of the theoretical uncertainty; however, we do 
not consider the result without fragmentation in our un- 
certainty estimate. 

Other perturbative QCD calculations are unlikely to 
give a much larger prompt neutrino flux than our upper 
limit, if saturation is indeed important. The theoretical 
expectation is that saturation is important at scales com- 
parable to jj, ~ m c for small x values. If it would turn 
out that saturation does not occur at the relevant energy 
scales, the flux is still not expected to be much larger than 
the PRS result. In Figure [SJ we therefore show a com- 
parison of the uncertainty (blue, dark band) compared to 
the proposed uncertainty range from Ref. [f| (magenta, 
light band) with their overlapping region (middle, light 
blue band). In this plot we take the NLO QCD result 
as the upper theoretical limit. This gives a larger up- 
per limit than in the earlier plots, which show only the 
uncertainty in the dipole model result. We stress that, 
since saturation is expected to be important on theoreti- 
cal grounds, the uncertainty band in Figure [5] is our main 
result. 

In order to obtain a flux as large as the upper line 
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log 10 E [GeV] 

FIG. 9: Prompt (solid line) and conventional (dashed lines) 
fluxes of v^+Vjj^ Ve + Ue, and u + +fJ>~ . The conventional fluxes 
are from Ref. fl5f | . The three prompt fluxes are approximately 
equal and are therefore here represented by the + flux. 



in the uncertainty band of Figure [8l we must multiply 
the upper uncertainty line of our DM result by a factor 
of 50. A cross section a factor of 50 times larger than 
the DM evaluation in proton-proton scattering would be 
incompatible with existing cross section measurements, 
as illustrated for example by Figure 6 of Ref. [3l[ , which 
compares the DM result for charm production to fixed- 
target experimental data. 

Measurable stau production rates from prompt atmo- 
spheric neutrinos as proposed in Ref. [6j would require 
the highest fluxes in the lighter band. Our evaluation of 
the prompt neutrino flux indicates that the upper limit 
of Ref. [fj is unrealistically large. The prompt neutrino 
flux is unlikely to be large enough for studying stau pro- 
duction from neutrino interactions with Earth and the 
subsequent detection in neutrino telescopes. 

The flavor decomposition of an atmospheric neutrino 
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FIG. 10: Prompt v r + v T flux (solid line) compared with the 
prompt + flux (dashed line). 
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FIG. 11: Dependence on zenith angle of prompt and conven- 
tional + fluxes. The solid and dashed lines are DM with 
fragmentation in the vertical and horizontal directions, re- 
spectively. The dotted and dash-dotted lines are the Gaisser- 
Honda conventional fluxes [ll| in the vertical and horizontal 
directions. 



tical conventional fluxes of muons, muon neutrinos and 
electron neutrinos (and their antiparticles) in Figure 
If experiments would be able to measure electron neu- 
trino fluxes, the prompt flux will start dominating over 
the conventional flux for much lower energy ~ 10 4 GeV 
than for the muon neutrino or muon fluxes. 

We note that the prompt flux of v T + v T from charm 
decays is much smaller than the other neutrino flavors 
[HI, since only the D s meson decays into v T . The v T + 
v r flux from D s decays is shown in Figure [TO] together 
with the prompt + flux. We have not included 
the contribution from i3-meson decays which could give 
a contribution on the order of 10-20% [l6| since i?-meson 
decays to v T plus tau are kinematically allowed. 

The vertical direction is the optimal direction for 
studying the prompt fluxes. In Figure [11] we show the 
prompt and conventional + fluxes in the vertical 
and horizontal directions. In the horizontal direction the 
prompt flux does not become larger than the conven- 
tional flux until very large energies ~ 10 7 GeV, where 
the actual number of neutrinos is quite low. 

In summary, we have computed prompt neutrino and 
muon fluxes from cosmic ray interactions in the atmo- 
sphere that produce charm pairs. Our evaluation of the 
fluxes takes parton saturation effects into account via the 
dipole model, a model with a parametric form guided by 
QCD and constrained by data. We find that saturation 
effects in the dipole model decrease the prompt fluxes 
above 10 5 GeV. Our estimate of the theoretical uncer- 
tainty in the predicted fluxes in the dipole model is on 
the order of a factor of two. In comparison to other QCD 
or dipole model evaluations of the prompt flux, the range 
of predictions is approximately a factor of 6. Future mea- 
surements of the high energy neutrino flux will provide 
interesting constraints on QCD-based evaluations of the 
prompt flux of neutrinos, however, the prompt neutrino 
flux is unlikely to be large enough to probe non-standard 
model interactions. 



signal may be an interesting way to explore the prompt 
contribution. The prompt neutrino fluxes of + Cy, 
and v e + z7 e are identical, since the charmed mesons de- 
cay equally likely into electrons or muons. The prompt 
// + + pT flux is approximately equal to the neutrino 
fluxes. However, this does not hold for the conven- 
tional fluxes. Charged pions decay almost exclusively 
into muons, so the muon neutrino and muon fluxes are 
much larger than the electron neutrino flux. We show the 
V[i + Up prompt flux together with the corresponding ver- 
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